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Abstract 


The time dependent governing acoustic- 
difference equations and boundary conditions are 
developed and solved for sound propagation In an 
axlsyametrlc (cylindrical) hard wall duct with a 
plug mean flow and spinning acoustic modes. The 
analysis begins with a harmonic sound tource radi- 
ating into a quieacent duct. This explicit iter- 
ation method then calculates stepwise In real time 
to obtain the transient as well as the "steady" 
state solutions of the acoustic field. The time 
dependent finite difference analysis has two advan- 
tages over the steady state finite difference and 
finite element techniques: (1) the elimination of 

large matrix storage requirements, and (2) shorter 
solution times under most conditions. 
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Nomenclature 

(ir 2 /^t)T ir M 

cell coefficient 

2(^r/^t) 2 r2n/C M 

cell coefficient 

ambient speed of sound, m/s 

coll coefficient sound 

cell coefficient 

(.5r 2 /^c5t)q r (l-M 2 )/C M 

cell coefficient 

frequency, Hz 

cell coefficient 

cell coefficient 

cell coefficient 

number of axial grid points 

/n 

cell coefficient 

number of transverse grid points 

Bessel function 

length of duct, m 

Mach number 

spinning mode number 

time-dependent dimensionless acoustic pressure, 
P(x,r, t) , P*/.*C* 2 

time-dependent dimensionless acoustic pressure 
associated with m mode 

time-dependent dimensionless acoustic pressure 



spatially dependent "steady" acoustic pressure, 
P(*,r) 

analytical solution for m spinning mode and 
n radial mode 


r radial coordinate, r*/r* 

r * radius of duct, m 

-3r radial grid spacing 

T* period, 1/f*, a 

t dimensionless time, t*/T* 

-5t time step 

l< axial acoustic velocity, u*/c* 

x axial coordinate, x*/rj 

axial grid spacing 
a eq. (18) 

°m n eigenvalue 
8 eq. (12) 

specific acoustic Impedance at exist 

Cm C e a+H/te> 

q r dimensionless frequency, r* f*/C* 

\ dimensionless axial wavelength 

t- angular coordinate 

,* ambient air density, kg/aP 

* eq. (19) 
angular frequency 

Subscripts 
e exit condition 

1 axial index (fig. 1) 

j radial index (fig. 1) 

k cell index 
m spinning mode nunfccr 

n radial mode number 

o ambient condition 
Superscripts 

* dimensional quantity 
k time step 

(1) real part 

(2) imaginary part 


Introduction 


Both finite difference and finite element 
numerical techniques have been developed to study 
sound propagation with axial variations in Mach 
number, wall impedance, and duct geometry. Gener- 
ally, the numerical solutions have been limited to 
low frequency sound and short ducts, because many 
grid points or elements were required to resolve 
the axial wavelength of the sound. For plane wave 
propagation, the number of axial grid points or 
'lements is directly proportional to the sound fre- 
quency and duct length* (eq. (77)), and inversely 
proportional to one minus the Mach number.* This 
later dependence also severely limits the nppli- 


1 


cation of numerical technique* for high Mach num- 
ber Inlet*. 

Custosurlly, the pre**ure and acoustic veloc- 
ities are aasuaed to be aiaple harmcnl function* 
of time. In this case, the time independent wave 
equation* can be uaed, 3 The aatricea associated 
with the numerical solution to the time Independe t 
equations (called steady state solutions herelr) 
must be solved exactly using such methods as Gauss 
ellmln Ion. As a result, large arrays of matrix 
elements must be stored. In an unpublished word 
at NASA Lewis Research Center by the author using 
reference 4, as well as In work of Qulmr (p. 3), 
the matrix has been modified to allow Iteration 
techniques; unfortunately, the convergence is too 
slow to be of any practical value. Other ap- 
proaches, such as in reference 6, might still 
allow the use of Iteration methods. 

Some special techniques have been developed to 
overcome the difficulties of high frequency, long 
ducts, and high Mach numbers . As shown In refer- 
ences 7 and 8, the wave envelope numerical technique 
can reduce the required number of grid points by an 
order of magnitude. In reference 9, this technique 
was used to optimize multielement liners of long 
lengths at high frequencies. At the present time, 
this technique has been applied only to the simple 
cases of no flow and plug flow. A numerical spatial 
marching technique was also developed in references 
10 and 11. Compared to the standard finite differ- 
ence or finite elesient boundary value approaches, 
the numerical marching technique Is orders of magni- 
tude shorter In computational cime and required 
computer storage. The marching technique is limited 
to high frequencies and to cases where reflections 
are small. 

A comprehensive literature suanary of other time 
independent finite difference and finite element 
techniques Is given elsewhere. 2 

As an alternative to the previously developed 
steady state theories, time dependent numerical so- 
lutions were developed for noise propagation in a 
two dimensional duct without flow* 2 and with par- 
allel sheared mean flow. ^ Advantageously, matrix 
storage requlrementa are significantly reduced In 
the time dependent snalysis. The analysis began 
with a noise source radiating into an Initially 
quiescent duct. Next, an explicit iteration 
method calculated stepwise in real time to ob- 
tain the transient as well as the "steady" state 
solution of the acoustic field. This time march- 
ing technique was found to be stable for both ro 
flow and plug flow. The time dependent analysis 
was found to be superior to the steady state finite 
difference and finite element techniques because of 
much shorter solution times and the reduction of 
matrix storage requirements. Also because matrix 
manipulation Is not required, the time-dependent 
approach is simpler to program and debug. 

In the present paper, the time dependent tech- 
nique developed In reference 12 will be extended 
to an axlsymnetric (cylindrical) geometry with 
spinning acoustic modes and plug flow. First, the 
governin'- acoustic equations and boundary conditions 
are presented. Next, the governing acoustic- 
difference equations are derived. Immediately 
following the mathematical development, numerical 
examples are presented and compared with the corre- 
sponding steady state analytical results. 


Governing Equations end Boundary Conditions 

The propagation of sound In an axiayamtetrlc 
cylindrical hard wall duct, as shown in figure 1, 
is described by the wave equation and appropriate 
impedance boundary conditions. 


Have Equations 


The wave equation In a circular duct with a 
mean flow can be expressed in dimensionless form as 
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d 2 ? 1 d 2 ? 1 1 dp 1 

dx 2 dr 2 r dr 


+ _L MLL . 2V , &L 

r 2 dd 2 ^ dxdt 


( 1 ) 


These and other symbols are defined in the nomen- 
clature. The dimensionless frequency q r is de- 
fined as 
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The asterisks denotes dimensional quantities. In 
many other references, the characteristic length 
used In the definition of q Is often the duct 
diameter. 

Because of the rotational nature of the rotor 
blades on a typical turbofan jet engine, large cir- 
cumferential variations In acoustic pressure will 
occur depending on blade number and engine rpm. A 
three dimensional solution for sound propagation, 
however, would be expensive to perform. Customar- 
ily, since the equations are linear, the circum- 
ferential acoustic pressure variations are decom- 
posed into spinning modes m: 

P’(*,r,6,t) - Ep.tx.r.tje 1 * (3) 

m 

The emanation Is over those modes that are likely 
to occur in a partlcal applications. Considering 
solutions with a single spinning lobe number m. 


P'(x,r,e,t) - P (x, r , t) e 1 ”^ (4) 


the wave equation (1) reduces to 



Equation (5) In difference form will be solved to 
determine the pressure in the duct. 


Hard Hall Boundary Condition 

The boundary condition at the surface of a 
hard wall duct la 


dp 

dr 
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( 6 ) 
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Entrance Condition 


The boundary condition at the source plane 
P(o,r,t) trill be eseuned to vary aa e lu * e * or In a 
dimenslonleaa fora aa e 12nt . Furthermore, the 
transverae preaaure variation (r-dlrectlon) will be 
aaauaed to correapond to the eigen functlone 
Jo(o m „r) aaaoclated with aode propagation In an 
Infinitely long hard wall duct. The elgenvaluea 
aaaoclated with mode (a,n) are tabulated In ref- 
erence 14, p. Sll and reference IS, p. 411. There- 
fore, the aource boundary condition uaed herein la 

P(°, r ,t) - J B (o„ „r)e i2nt (7) 


Finally, the ten d^P/dadt in Che wave 
equation will alao be evaluated at the duct exit. 
Differentiating equation (10) with tlam yielda: 
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(13) 


Centerline Condition 


The difference equation will be developed for 
taodea with m equal or greater than 1. In these 
caaea, the preaaure at the centerline la equal to 
aero: 


Exit Impedance 


P(x,o,t) - 0 m > 1 (14) 


The boundary condition at the exit of the duct 
can be expreaaed In terms of a specific acoustic 
impedance defined as 


■ u 


( 8 ) 


The equation governing the acoustic velocity U la 
th usual linearized momentum equation: 


du 

dt 


1 dp 
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M_ du 
r 


(9) 


Substituting equation (8) Into equation (9) yielda: 
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Initial Condition 

For tlmea equal to or leas than zero, the duct 
la assumed quiescent, that la, the acoustic pressure 
la taken to be zero. For times greater than zero, 
the application of Che noise aource (eq. (7)) will 
drive the pressures In the duct. 


Complex Notation 

Because of the introduction of complex nota- 
tion for the noise source, the acoustic pressure P 
is complex. The superscript (1) will represent the 
real term while (2) will represent the imaginary 
term. 


P - p(l) + IP (2) 


(15) 


In equation (10), C e has been assumed to be the 
impedance associated with mode propagation down an 
infinite duct. For transmission ot a single acous- 
tic mode without reflection, the exit Impedance 
i S 16 


1 - MB 
p - M 


(ID 


where 


Difference Equations 

Instead of a continuous solution In space and 
time, the finite-difference approximations will 
determine the pressure at Isolated grid points In 
space as shown in figure 1, and at discrete time 
steps At. Starting from the known Initial condi- 
tions at t - 0 and the boundary conditions , the 
finite-difference algorithm will march-out the solu- 
tion to later times. 



Centrsl Region (Cell No. 1) 

Away from the duct boundaries. In Cell no. 1 
of figure 1, the second and first derivatives In the 
wave equation (eq. (5)) can be represented by the 
usual central differences In time and spaced” 

(p. 884) 
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where 1 and j denote the space Indices, k the 
tlite laden, aad Ac, Ac, aad At are the space and 
tine aesh spacing, respectively, the tlmejasao- 
clated with equation (16) Is defined as t** 1 
- r + it ■ (k + l)At. Equation (16) can be 
written in the form 

*?.] * <j + T t\ ■ «- 1} < i7 > 


Boundary Condition (Cells Ho. 2 to 6) 

The erpreaslons for the difference equations 
at the boundaries are conpllcated by the impedance 
condition and the change In area of cells no. 2 to 
6 In figure 1. The governing difference equations 
can be developed by an Integration process In which 
the wave equation (eq. (5)) Is Integrated over the 
area of the cells and time: 


where 



( 18 ) 


/ t+At/2 a, 

li 

At/2 •'"'Cell 


dxdrdt 


In this case, T. represents all the remaining dlf- 
ference terms on the right hand side of the equa- 
tion. The subscript q stands for Che cells shown 
In figure 1. In this case, q equals one for the 
central cell. 

For ease In bookkeeping, I Is broken down as 
follows 



area 


rf) + + 1*1 

dx? clr2 r dr 


% ■ Vi-i.J + Vi.j-1 + v*.j + V*.j+i 


+ e q P k+l.J + f q P ltj 


* ♦ vn.i * vw., 


(19) 


Substituting equation (19) Into equation (17) and 
solving for the pressure plf'j, yields 

'“i ■ * " A) -‘ 

_k . _k , „k | 

* P i,J Vi,J+i e q p i+i,JJ 

* [vEi.j + ViJi.J < 20) 


The coefficients are listed In tsble I. 

Equation (20) Is an algorithm which permits 
marching out solutions from known values of pressure 
at times associated with k and k-1. The term 
Pjtj i is also a known quantity. At 1 * 2, the 
valueof Is known from the source condition, 

equation (7). At 1-3, the value of Is 

known from the previous calculation at 1 * J 2. 

Thus, the solution starts at the spatial position 
1-2 and proceeda outward to 1 - I. The proce- 
dure Is explicit since ell the past values of P k 
are known as the new values of k + 1 are computed. 
For the special case at t - 0, the values of the 
pressure associated with the k-1 value are zero 
from the assumed Initial condition. 


- =7 P - 2.)M — Mdxdrdt (21) 

r* oxo tj 

The procedure for the temporal end spatial Integra- 
tion over the cell area is documented in refer- 
ence 12, and therefore will not be presented herein. 
However, some extra steps that were not presented 
In reference 12 are listed in the Appendix of this 
paper. 

The finite-difference approximation for the 
various cells shown In figure 1 are identical to 
equation (20) with the various cell coefficients 
listed in table I. Cells 3 and 6 have not been 
evaluated since the pressure at the centerline is 
equal to zero for spinning modes . 


Spatial Mesh Size 

The mesh spacing Ax and Ar must be re- 
stricted to small values to reduce the truncation 
error. To resolve the oscillatory nature of the 
pressure, the required number of grid points in the 
axial direction suggested In reference 1 was 

12 L*/r* 

I > l (22) 

where A Is the axial wavelength of the sound. For 
spinning mode propagation In a seml-inftnlte hard 
cylindrical duct, A equals 

1 - 

A - — (23) 

T) r (p - M) 

therefore 

I > 12i]_L*/r* & ' (24) 

(1 - M 2 ) 


Similarly. Che number of points used In the 
r direction is*2 


J - 12, ir 


(25) 
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Stability 

In the explicit time marching approach used 
here, round-off error* can grow In an unbounded 
fashion and destroy the solution If the time Incre- 
ment .it Is too large or the iteration scheme is 
Improperly posed. In the present paper, numerical 
experimentation was used to determine the time In- 
crement .it for which the explicit solution re- 
mained stable. The Initial guess for the time in- 
crement was baaed on the following equation devel- 
oped In reference 13 for plane wave propagation In 
straight ducts: 


•'* < ’> r 


- ImI) 
Vl + Oir/Ax) 2 


(26) 


Steady State Pressures 


In the sample problems to be presented in the 
next section, the time dependent results will be 
compared to the results of the steady harmonic 
solutions (p. 509). The puipusO of this section 
is to show tiie rational for constructing a steady 
state solution from the time-dependent results. 


Steady Harmonic Solution 

The steady acoustic pressure p(x,r) is defined 
as the solution to equation (5) when the pressure 
is assumed to be a simple harmonic function of time. 
For a semi-infinite hard wall circular duct with 
transmission of a single spinning mode, the analyti- 
cal solution for p (x,r) is 
m n 


of the steady harmonic solution given by equa- 
tion (27). 


Zero Mach Number 

Numerical and analytical values of the "steady" 
pressure p(x,r) are computed for the case of a 
hard wall infinite duct with an m • 3 spinning 
mode. The calculation was made with a length to 
radius ratio of 1 and dimensionless frequencies of 
q r * 2, 1, and 0.7. The analytical and numerical 
values of the acoustic pressure profiles along the 
duct wall (r ' l) are shown in figures 2(a), (b), 
and (c). As seen in figure 2, the steady and time- 
dependent analyses are in good agreement. 

The cut-off frequency for the (3,0) mode is 
0.66854. For <j r » 0.5, the numerical solution as 
formulated did not converge to the analytical solu- 
tion. The numerical solution oscillated about the 
analytical solution. Therefore, at the present 
time, the numerical procedure can only be applied 
to propagating acoustic modes. 


Finite Mach Number 


As another example of the time-dependent anal- 
ysis, numerical and analytical values of the 
"steady" pressure p(x,r) are computed at « 1 
but with uniform Mach numb./s of +0.5 and -0.5. 

The analytical and numerical values o.~ the acoustic 
pressure are shown in figures 3(a) and (b). As 
seen in figure 3, agreement between Che analytical 
and numerical theory is good. 


Conclusions 


-i2m, r (;'--M)x 

Pm n< x - r > = J m (a s> n r > c 1_m2 (* > °> 

(27) 

where g is given by equation (12). In the next 
section of the present paper, the transient solution 
to the semi-infinite hard wall duct will be com- 
pared to equation (27). 


• anslent Solution 

Although multiple values of pressure are 
calculated at each time step, only the latest value 
values (V, k-1) need to be stored. After the ini- 
tial transient has died out (checked numerically) , 
the time dependent results can be compared to the 
"steady" state results (eq. (27)) simply by divid- 
ing by e^ nt , that is, 


A time dependent numerical procedure was de- 
veloped for no flow and plug flow in a cylindrical 
duct. This explicit time marching technique was 
found to be stable for propagating acoustic modes 
with or without flow. For nonpropagating (cut-off) 
acoustic modes, however, the time dependent solu- 
tion oscillates about the correct solution. This 
will ultimately have to be resolved since cut-off 
modes vill be encountered in real problems. In the 
mean time, the steady state numerical techniques 
can be conveniently employed near cut-off. Be- 
cause of the relatively long axial wave length near 
cut-off (see fig. 2(c)), the required number of 
grid points (or elements) is relatively small. 

By eliminating large matrix storage require- 
ments, numerical calculations of high sound fre- 
quencies in turbojet inlets are now possible. Be- 
cause manipulation of matrices is omitted, the time 
dependent approach is relatively easy to program 
and debug. 


p(x,r) - (28) 

A2nt 
e 

The bookkeeping and graphical output is held to a 
minimum by the use of equation (28). 

Sample Calculations 

In two sample problems to follow, the time- 
dependent results will be compared to the results 


Appendix 

Finite-Difference Equations at Boundary 

The derivation of the difference equations for 
the boundary cells was presented in reference 12. 
However, some additional terms appear in this paper 
which were not covered in reference 12. The treat- 
ment of these new terms are now considered. 


5 




References 

1. Baumelster, K. J. and Bittner, E. C. , "Numerical 

Simulation of Noise Piopagatlon In Jet En- 
gine Ducts," NASA TN D-7339, 1973. 

2. Baumelster, K. J. and Rice, E. J., "A Difference 

Theory for Noise Propagation In an Acoustical- 
ly Lined Duct with Mean Flow," A eroacousttcs : 
Jet and Combustion Noise: Duct Acoustics , 
edited by H. T. Nagamatsu, J. V. O'Keefe, and 
I. R. Schwartz, Progress in Astronautics and 
Aeronautics Series, Vol. 37, American Insti- 
tute of Aeronautics and Astronautics, New 
York, 1975, pp. 435-953. 

3. Goldstein, M. E., Aeroacoustlcs . McGraw-Hill, 

New York, 1976. 

4. Beaublen, M. J. and Uexler, A., "Iterative, 

Finite Difference Solution of Interior Eigen- 
values and Eigenfunctions of Laplace's Opera- 
tor," Corajniter_jnurnal, Vol. 14, Aug, 1971, 
pp. 263-269. 

5. Quinn, D. W., "A Finite Element Method for 

Computing Sound Propagation in Ducts Contain- 
ing Flow," AIAA Paper 79-0661, Mar. 1979. 

6. Browne, B. T. and Lawrenson, P. J., "Numerical 

Solution of an Elliptic Boundary-Value Problem 
In the Complex Variable," Institute of Mathe- 
matics and its Applications, Vol. 17, (3), 
1976, pp. 311-327. 

7. Baumelster, K. J., "Analysis of Sound Propaga- 

tion in Ducts Using the Wave Envelope Con- 
cept," NASA TN D-7719, 1974. 

8. Baumelster, K. J., "Finite-Difference Theory for 

Sound Propagation in a Lined Duct with Uniform 
Flow Using the Wave Envelope Concept," NASA 
TP- 1001, 1977. 

9. Baumelster, K. J., "Optimized Multisectioned 

Acoustic Liners," AIAA Paper 79-0182, Jan. 
1979. 


10. Baumelster, K. J., "Numerical Spatial Marching 

Techniques for Estimating Duct Attenuation 
and Source Pressure Profiles," NASA TM-78857, 

1978. 

11. Baumelster, K. J., "Numerical Spatial Marching 

Techniques in Duct Acoustics," Acoustical 
Society of America, Journal, Vol. 65, Feb. 

1979, pp. 297-306. 

12. Baumelster, K. J., "Time Dependent Difference 

Theory for Noise Propagation in Jet Engine 
Ducts," AIAA Paper 80-0098, Jan. 1980. 

13. Baumelster, K. J., "A Time Deperdent Difference 

Theory for Sound Propagation in Ducts with 
Flow," NASA TM-79302 , 1979. 

14. Morse, P. M. and Ingard, K. U. , Theoretical 

Acoustics . McGraw Hill, New York, 1968. 

15. Abrsmwltz, M. and Stegun, I. A., eds., Handbook 

of Mathematical Functions with Formulas , 
Graphs and Mathematical Tables . National 
Bureau of Standards, Applied Mathematics 
Series, 55, 1964. 

16. Tag, I. A. and Lumsdalne, E., "An efficient 

Finite Element Technique for Sound Propaga- 
tion in Axisymnetrlc hard Wall Ducts Carrying 
High Subsonic Mach Number Flows," AIAA Paper 
78-1154, July 1978. 


6 


TABLE I. - CELL SUMMARY 




Si 

PLANE -n 





DIMENSIONLESS ACOUSTIC PRESSURE, p 




lb) FREQUENCY n r ■ L 0, At ■ a 00444, t ■ 3. 006. 



(0 FREQUENCY r? r ■ Q, 7, At ■ a 003108, t ■ 10. 7a 

Figure 2. - Analytical and numerical pressure profiles 
for spinning wave propagation in an infinite hard wall 
duct (m ■ 3, n ■ 0, a K -4.20119, J ■ 20, r ■ 1). 



DIMENSIONLESS ACOUSTIC PRESSURE, p 




(b) MACH NUMBER M ■ -0. 5, £t m 0. 001631, t - 4.002. 

Figure 3. - Analytical and numerical pressure profiles 
for spinning wave propagation in an infinite hard wall 
duct with flow (m « 3, n ■ 0, 030 ■ 4.20119, J • 20, 
r • 1 , and 7 j r ■ 1 ). 


